// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPARSEMATRIX_H
#define EIGEN_SPARSEMATRIX_H

namespace Eigen {

/** \ingroup SparseCore_Module
  *
  * \class SparseMatrix
  *
  * \brief A versatible sparse matrix representation
  *
  * This class implements a more versatile variants of the common \em compressed row/column storage format.
  * Each colmun's (resp. row) non zeros are stored as a pair of value with associated row (resp. colmiun) index.
  * All the non zeros are stored in a single large buffer. Unlike the \em compressed format, there might be extra
  * space in between the nonzeros of two successive colmuns (resp. rows) such that insertion of new non-zero
  * can be done with limited memory reallocation and copies.
  *
  * A call to the function makeCompressed() turns the matrix into the standard \em compressed format
  * compatible with many library.
  *
  * More details on this storage sceheme are given in the \ref TutorialSparse "manual pages".
  *
  * \tparam _Scalar the scalar type, i.e. the type of the coefficients
  * \tparam _Options Union of bit flags controlling the storage scheme. Currently the only possibility
  *                 is ColMajor or RowMajor. The default is 0 which means column-major.
  * \tparam _StorageIndex the type of the indices. It has to be a \b signed type (e.g., short, int, std::ptrdiff_t). Default is \c int.
  *
  * \warning In %Eigen 3.2, the undocumented type \c SparseMatrix::Index was improperly defined as the storage index type (e.g., int),
  *          whereas it is now (starting from %Eigen 3.3) deprecated and always defined as Eigen::Index.
  *          Codes making use of \c SparseMatrix::Index, might thus likely have to be changed to use \c SparseMatrix::StorageIndex instead.
  *
  * This class can be extended with the help of the plugin mechanism described on the page
  * \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_SPARSEMATRIX_PLUGIN.
  */

namespace internal {
    template <typename _Scalar, int _Options, typename _StorageIndex> struct traits<SparseMatrix<_Scalar, _Options, _StorageIndex>>
    {
        typedef _Scalar Scalar;
        typedef _StorageIndex StorageIndex;
        typedef Sparse StorageKind;
        typedef MatrixXpr XprKind;
        enum
        {
            RowsAtCompileTime = Dynamic,
            ColsAtCompileTime = Dynamic,
            MaxRowsAtCompileTime = Dynamic,
            MaxColsAtCompileTime = Dynamic,
            Flags = _Options | NestByRefBit | LvalueBit | CompressedAccessBit,
            SupportedAccessPatterns = InnerRandomAccessPattern
        };
    };

    template <typename _Scalar, int _Options, typename _StorageIndex, int DiagIndex>
    struct traits<Diagonal<SparseMatrix<_Scalar, _Options, _StorageIndex>, DiagIndex>>
    {
        typedef SparseMatrix<_Scalar, _Options, _StorageIndex> MatrixType;
        typedef typename ref_selector<MatrixType>::type MatrixTypeNested;
        typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested;

        typedef _Scalar Scalar;
        typedef Dense StorageKind;
        typedef _StorageIndex StorageIndex;
        typedef MatrixXpr XprKind;

        enum
        {
            RowsAtCompileTime = Dynamic,
            ColsAtCompileTime = 1,
            MaxRowsAtCompileTime = Dynamic,
            MaxColsAtCompileTime = 1,
            Flags = LvalueBit
        };
    };

    template <typename _Scalar, int _Options, typename _StorageIndex, int DiagIndex>
    struct traits<Diagonal<const SparseMatrix<_Scalar, _Options, _StorageIndex>, DiagIndex>>
        : public traits<Diagonal<SparseMatrix<_Scalar, _Options, _StorageIndex>, DiagIndex>>
    {
        enum
        {
            Flags = 0
        };
    };

}  // end namespace internal

template <typename _Scalar, int _Options, typename _StorageIndex>
class SparseMatrix : public SparseCompressedBase<SparseMatrix<_Scalar, _Options, _StorageIndex>>
{
    typedef SparseCompressedBase<SparseMatrix> Base;
    using Base::convert_index;
    friend class SparseVector<_Scalar, 0, _StorageIndex>;
    template <typename, typename, typename, typename, typename> friend struct internal::Assignment;

public:
    using Base::isCompressed;
    using Base::nonZeros;
    EIGEN_SPARSE_PUBLIC_INTERFACE(SparseMatrix)
    using Base::operator+=;
    using Base::operator-=;

    typedef MappedSparseMatrix<Scalar, Flags> Map;
    typedef Diagonal<SparseMatrix> DiagonalReturnType;
    typedef Diagonal<const SparseMatrix> ConstDiagonalReturnType;
    typedef typename Base::InnerIterator InnerIterator;
    typedef typename Base::ReverseInnerIterator ReverseInnerIterator;

    using Base::IsRowMajor;
    typedef internal::CompressedStorage<Scalar, StorageIndex> Storage;
    enum
    {
        Options = _Options
    };

    typedef typename Base::IndexVector IndexVector;
    typedef typename Base::ScalarVector ScalarVector;

protected:
    typedef SparseMatrix<Scalar, (Flags & ~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0)> TransposedSparseMatrix;

    Index m_outerSize;
    Index m_innerSize;
    StorageIndex* m_outerIndex;
    StorageIndex* m_innerNonZeros;  // optional, if null then the data is compressed
    Storage m_data;

public:
    /** \returns the number of rows of the matrix */
    inline Index rows() const { return IsRowMajor ? m_outerSize : m_innerSize; }
    /** \returns the number of columns of the matrix */
    inline Index cols() const { return IsRowMajor ? m_innerSize : m_outerSize; }

    /** \returns the number of rows (resp. columns) of the matrix if the storage order column major (resp. row major) */
    inline Index innerSize() const { return m_innerSize; }
    /** \returns the number of columns (resp. rows) of the matrix if the storage order column major (resp. row major) */
    inline Index outerSize() const { return m_outerSize; }

    /** \returns a const pointer to the array of values.
      * This function is aimed at interoperability with other libraries.
      * \sa innerIndexPtr(), outerIndexPtr() */
    inline const Scalar* valuePtr() const { return m_data.valuePtr(); }
    /** \returns a non-const pointer to the array of values.
      * This function is aimed at interoperability with other libraries.
      * \sa innerIndexPtr(), outerIndexPtr() */
    inline Scalar* valuePtr() { return m_data.valuePtr(); }

    /** \returns a const pointer to the array of inner indices.
      * This function is aimed at interoperability with other libraries.
      * \sa valuePtr(), outerIndexPtr() */
    inline const StorageIndex* innerIndexPtr() const { return m_data.indexPtr(); }
    /** \returns a non-const pointer to the array of inner indices.
      * This function is aimed at interoperability with other libraries.
      * \sa valuePtr(), outerIndexPtr() */
    inline StorageIndex* innerIndexPtr() { return m_data.indexPtr(); }

    /** \returns a const pointer to the array of the starting positions of the inner vectors.
      * This function is aimed at interoperability with other libraries.
      * \sa valuePtr(), innerIndexPtr() */
    inline const StorageIndex* outerIndexPtr() const { return m_outerIndex; }
    /** \returns a non-const pointer to the array of the starting positions of the inner vectors.
      * This function is aimed at interoperability with other libraries.
      * \sa valuePtr(), innerIndexPtr() */
    inline StorageIndex* outerIndexPtr() { return m_outerIndex; }

    /** \returns a const pointer to the array of the number of non zeros of the inner vectors.
      * This function is aimed at interoperability with other libraries.
      * \warning it returns the null pointer 0 in compressed mode */
    inline const StorageIndex* innerNonZeroPtr() const { return m_innerNonZeros; }
    /** \returns a non-const pointer to the array of the number of non zeros of the inner vectors.
      * This function is aimed at interoperability with other libraries.
      * \warning it returns the null pointer 0 in compressed mode */
    inline StorageIndex* innerNonZeroPtr() { return m_innerNonZeros; }

    /** \internal */
    inline Storage& data() { return m_data; }
    /** \internal */
    inline const Storage& data() const { return m_data; }

    /** \returns the value of the matrix at position \a i, \a j
      * This function returns Scalar(0) if the element is an explicit \em zero */
    inline Scalar coeff(Index row, Index col) const
    {
        eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols());

        const Index outer = IsRowMajor ? row : col;
        const Index inner = IsRowMajor ? col : row;
        Index end = m_innerNonZeros ? m_outerIndex[outer] + m_innerNonZeros[outer] : m_outerIndex[outer + 1];
        return m_data.atInRange(m_outerIndex[outer], end, StorageIndex(inner));
    }

    /** \returns a non-const reference to the value of the matrix at position \a i, \a j
      *
      * If the element does not exist then it is inserted via the insert(Index,Index) function
      * which itself turns the matrix into a non compressed form if that was not the case.
      *
      * This is a O(log(nnz_j)) operation (binary search) plus the cost of insert(Index,Index)
      * function if the element does not already exist.
      */
    inline Scalar& coeffRef(Index row, Index col)
    {
        eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols());

        const Index outer = IsRowMajor ? row : col;
        const Index inner = IsRowMajor ? col : row;

        Index start = m_outerIndex[outer];
        Index end = m_innerNonZeros ? m_outerIndex[outer] + m_innerNonZeros[outer] : m_outerIndex[outer + 1];
        eigen_assert(end >= start && "you probably called coeffRef on a non finalized matrix");
        if (end <= start)
            return insert(row, col);
        const Index p = m_data.searchLowerIndex(start, end - 1, StorageIndex(inner));
        if ((p < end) && (m_data.index(p) == inner))
            return m_data.value(p);
        else
            return insert(row, col);
    }

    /** \returns a reference to a novel non zero coefficient with coordinates \a row x \a col.
      * The non zero coefficient must \b not already exist.
      *
      * If the matrix \c *this is in compressed mode, then \c *this is turned into uncompressed
      * mode while reserving room for 2 x this->innerSize() non zeros if reserve(Index) has not been called earlier.
      * In this case, the insertion procedure is optimized for a \e sequential insertion mode where elements are assumed to be
      * inserted by increasing outer-indices.
      * 
      * If that's not the case, then it is strongly recommended to either use a triplet-list to assemble the matrix, or to first
      * call reserve(const SizesType &) to reserve the appropriate number of non-zero elements per inner vector.
      *
      * Assuming memory has been appropriately reserved, this function performs a sorted insertion in O(1)
      * if the elements of each inner vector are inserted in increasing inner index order, and in O(nnz_j) for a random insertion.
      *
      */
    Scalar& insert(Index row, Index col);

public:
    /** Removes all non zeros but keep allocated memory
      *
      * This function does not free the currently allocated memory. To release as much as memory as possible,
      * call \code mat.data().squeeze(); \endcode after resizing it.
      * 
      * \sa resize(Index,Index), data()
      */
    inline void setZero()
    {
        m_data.clear();
        memset(m_outerIndex, 0, (m_outerSize + 1) * sizeof(StorageIndex));
        if (m_innerNonZeros)
            memset(m_innerNonZeros, 0, (m_outerSize) * sizeof(StorageIndex));
    }

    /** Preallocates \a reserveSize non zeros.
      *
      * Precondition: the matrix must be in compressed mode. */
    inline void reserve(Index reserveSize)
    {
        eigen_assert(isCompressed() && "This function does not make sense in non compressed mode.");
        m_data.reserve(reserveSize);
    }

#ifdef EIGEN_PARSED_BY_DOXYGEN
    /** Preallocates \a reserveSize[\c j] non zeros for each column (resp. row) \c j.
      *
      * This function turns the matrix in non-compressed mode.
      * 
      * The type \c SizesType must expose the following interface:
        \code
        typedef value_type;
        const value_type& operator[](i) const;
        \endcode
      * for \c i in the [0,this->outerSize()[ range.
      * Typical choices include std::vector<int>, Eigen::VectorXi, Eigen::VectorXi::Constant, etc.
      */
    template <class SizesType> inline void reserve(const SizesType& reserveSizes);
#else
    template <class SizesType>
    inline void reserve(const SizesType& reserveSizes,
                        const typename SizesType::value_type& enableif =
#if (!EIGEN_COMP_MSVC) || (EIGEN_COMP_MSVC >= 1500)  // MSVC 2005 fails to compile with this typename
                            typename
#endif
                        SizesType::value_type())
    {
        EIGEN_UNUSED_VARIABLE(enableif);
        reserveInnerVectors(reserveSizes);
    }
#endif  // EIGEN_PARSED_BY_DOXYGEN
protected:
    template <class SizesType> inline void reserveInnerVectors(const SizesType& reserveSizes)
    {
        if (isCompressed())
        {
            Index totalReserveSize = 0;
            // turn the matrix into non-compressed mode
            m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
            if (!m_innerNonZeros)
                internal::throw_std_bad_alloc();

            // temporarily use m_innerSizes to hold the new starting points.
            StorageIndex* newOuterIndex = m_innerNonZeros;

            StorageIndex count = 0;
            for (Index j = 0; j < m_outerSize; ++j)
            {
                newOuterIndex[j] = count;
                count += reserveSizes[j] + (m_outerIndex[j + 1] - m_outerIndex[j]);
                totalReserveSize += reserveSizes[j];
            }
            m_data.reserve(totalReserveSize);
            StorageIndex previousOuterIndex = m_outerIndex[m_outerSize];
            for (Index j = m_outerSize - 1; j >= 0; --j)
            {
                StorageIndex innerNNZ = previousOuterIndex - m_outerIndex[j];
                for (Index i = innerNNZ - 1; i >= 0; --i)
                {
                    m_data.index(newOuterIndex[j] + i) = m_data.index(m_outerIndex[j] + i);
                    m_data.value(newOuterIndex[j] + i) = m_data.value(m_outerIndex[j] + i);
                }
                previousOuterIndex = m_outerIndex[j];
                m_outerIndex[j] = newOuterIndex[j];
                m_innerNonZeros[j] = innerNNZ;
            }
            if (m_outerSize > 0)
                m_outerIndex[m_outerSize] = m_outerIndex[m_outerSize - 1] + m_innerNonZeros[m_outerSize - 1] + reserveSizes[m_outerSize - 1];

            m_data.resize(m_outerIndex[m_outerSize]);
        }
        else
        {
            StorageIndex* newOuterIndex = static_cast<StorageIndex*>(std::malloc((m_outerSize + 1) * sizeof(StorageIndex)));
            if (!newOuterIndex)
                internal::throw_std_bad_alloc();

            StorageIndex count = 0;
            for (Index j = 0; j < m_outerSize; ++j)
            {
                newOuterIndex[j] = count;
                StorageIndex alreadyReserved = (m_outerIndex[j + 1] - m_outerIndex[j]) - m_innerNonZeros[j];
                StorageIndex toReserve = std::max<StorageIndex>(reserveSizes[j], alreadyReserved);
                count += toReserve + m_innerNonZeros[j];
            }
            newOuterIndex[m_outerSize] = count;

            m_data.resize(count);
            for (Index j = m_outerSize - 1; j >= 0; --j)
            {
                Index offset = newOuterIndex[j] - m_outerIndex[j];
                if (offset > 0)
                {
                    StorageIndex innerNNZ = m_innerNonZeros[j];
                    for (Index i = innerNNZ - 1; i >= 0; --i)
                    {
                        m_data.index(newOuterIndex[j] + i) = m_data.index(m_outerIndex[j] + i);
                        m_data.value(newOuterIndex[j] + i) = m_data.value(m_outerIndex[j] + i);
                    }
                }
            }

            std::swap(m_outerIndex, newOuterIndex);
            std::free(newOuterIndex);
        }
    }

public:
    //--- low level purely coherent filling ---

    /** \internal
      * \returns a reference to the non zero coefficient at position \a row, \a col assuming that:
      * - the nonzero does not already exist
      * - the new coefficient is the last one according to the storage order
      *
      * Before filling a given inner vector you must call the statVec(Index) function.
      *
      * After an insertion session, you should call the finalize() function.
      *
      * \sa insert, insertBackByOuterInner, startVec */
    inline Scalar& insertBack(Index row, Index col) { return insertBackByOuterInner(IsRowMajor ? row : col, IsRowMajor ? col : row); }

    /** \internal
      * \sa insertBack, startVec */
    inline Scalar& insertBackByOuterInner(Index outer, Index inner)
    {
        eigen_assert(Index(m_outerIndex[outer + 1]) == m_data.size() && "Invalid ordered insertion (invalid outer index)");
        eigen_assert((m_outerIndex[outer + 1] - m_outerIndex[outer] == 0 || m_data.index(m_data.size() - 1) < inner) &&
                     "Invalid ordered insertion (invalid inner index)");
        Index p = m_outerIndex[outer + 1];
        ++m_outerIndex[outer + 1];
        m_data.append(Scalar(0), inner);
        return m_data.value(p);
    }

    /** \internal
      * \warning use it only if you know what you are doing */
    inline Scalar& insertBackByOuterInnerUnordered(Index outer, Index inner)
    {
        Index p = m_outerIndex[outer + 1];
        ++m_outerIndex[outer + 1];
        m_data.append(Scalar(0), inner);
        return m_data.value(p);
    }

    /** \internal
      * \sa insertBack, insertBackByOuterInner */
    inline void startVec(Index outer)
    {
        eigen_assert(m_outerIndex[outer] == Index(m_data.size()) && "You must call startVec for each inner vector sequentially");
        eigen_assert(m_outerIndex[outer + 1] == 0 && "You must call startVec for each inner vector sequentially");
        m_outerIndex[outer + 1] = m_outerIndex[outer];
    }

    /** \internal
      * Must be called after inserting a set of non zero entries using the low level compressed API.
      */
    inline void finalize()
    {
        if (isCompressed())
        {
            StorageIndex size = internal::convert_index<StorageIndex>(m_data.size());
            Index i = m_outerSize;
            // find the last filled column
            while (i >= 0 && m_outerIndex[i] == 0) --i;
            ++i;
            while (i <= m_outerSize)
            {
                m_outerIndex[i] = size;
                ++i;
            }
        }
    }

    //---

    template <typename InputIterators> void setFromTriplets(const InputIterators& begin, const InputIterators& end);

    template <typename InputIterators, typename DupFunctor> void setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func);

    void sumupDuplicates() { collapseDuplicates(internal::scalar_sum_op<Scalar, Scalar>()); }

    template <typename DupFunctor> void collapseDuplicates(DupFunctor dup_func = DupFunctor());

    //---

    /** \internal
      * same as insert(Index,Index) except that the indices are given relative to the storage order */
    Scalar& insertByOuterInner(Index j, Index i) { return insert(IsRowMajor ? j : i, IsRowMajor ? i : j); }

    /** Turns the matrix into the \em compressed format.
      */
    void makeCompressed()
    {
        if (isCompressed())
            return;

        eigen_internal_assert(m_outerIndex != 0 && m_outerSize > 0);

        Index oldStart = m_outerIndex[1];
        m_outerIndex[1] = m_innerNonZeros[0];
        for (Index j = 1; j < m_outerSize; ++j)
        {
            Index nextOldStart = m_outerIndex[j + 1];
            Index offset = oldStart - m_outerIndex[j];
            if (offset > 0)
            {
                for (Index k = 0; k < m_innerNonZeros[j]; ++k)
                {
                    m_data.index(m_outerIndex[j] + k) = m_data.index(oldStart + k);
                    m_data.value(m_outerIndex[j] + k) = m_data.value(oldStart + k);
                }
            }
            m_outerIndex[j + 1] = m_outerIndex[j] + m_innerNonZeros[j];
            oldStart = nextOldStart;
        }
        std::free(m_innerNonZeros);
        m_innerNonZeros = 0;
        m_data.resize(m_outerIndex[m_outerSize]);
        m_data.squeeze();
    }

    /** Turns the matrix into the uncompressed mode */
    void uncompress()
    {
        if (m_innerNonZeros != 0)
            return;
        m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
        for (Index i = 0; i < m_outerSize; i++) { m_innerNonZeros[i] = m_outerIndex[i + 1] - m_outerIndex[i]; }
    }

    /** Suppresses all nonzeros which are \b much \b smaller \b than \a reference under the tolerance \a epsilon */
    void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits<RealScalar>::dummy_precision())
    {
        prune(default_prunning_func(reference, epsilon));
    }

    /** Turns the matrix into compressed format, and suppresses all nonzeros which do not satisfy the predicate \a keep.
      * The functor type \a KeepFunc must implement the following function:
      * \code
      * bool operator() (const Index& row, const Index& col, const Scalar& value) const;
      * \endcode
      * \sa prune(Scalar,RealScalar)
      */
    template <typename KeepFunc> void prune(const KeepFunc& keep = KeepFunc())
    {
        // TODO optimize the uncompressed mode to avoid moving and allocating the data twice
        makeCompressed();

        StorageIndex k = 0;
        for (Index j = 0; j < m_outerSize; ++j)
        {
            Index previousStart = m_outerIndex[j];
            m_outerIndex[j] = k;
            Index end = m_outerIndex[j + 1];
            for (Index i = previousStart; i < end; ++i)
            {
                if (keep(IsRowMajor ? j : m_data.index(i), IsRowMajor ? m_data.index(i) : j, m_data.value(i)))
                {
                    m_data.value(k) = m_data.value(i);
                    m_data.index(k) = m_data.index(i);
                    ++k;
                }
            }
        }
        m_outerIndex[m_outerSize] = k;
        m_data.resize(k, 0);
    }

    /** Resizes the matrix to a \a rows x \a cols matrix leaving old values untouched.
      *
      * If the sizes of the matrix are decreased, then the matrix is turned to \b uncompressed-mode
      * and the storage of the out of bounds coefficients is kept and reserved.
      * Call makeCompressed() to pack the entries and squeeze extra memory.
      *
      * \sa reserve(), setZero(), makeCompressed()
      */
    void conservativeResize(Index rows, Index cols)
    {
        // No change
        if (this->rows() == rows && this->cols() == cols)
            return;

        // If one dimension is null, then there is nothing to be preserved
        if (rows == 0 || cols == 0)
            return resize(rows, cols);

        Index innerChange = IsRowMajor ? cols - this->cols() : rows - this->rows();
        Index outerChange = IsRowMajor ? rows - this->rows() : cols - this->cols();
        StorageIndex newInnerSize = convert_index(IsRowMajor ? cols : rows);

        // Deals with inner non zeros
        if (m_innerNonZeros)
        {
            // Resize m_innerNonZeros
            StorageIndex* newInnerNonZeros = static_cast<StorageIndex*>(std::realloc(m_innerNonZeros, (m_outerSize + outerChange) * sizeof(StorageIndex)));
            if (!newInnerNonZeros)
                internal::throw_std_bad_alloc();
            m_innerNonZeros = newInnerNonZeros;

            for (Index i = m_outerSize; i < m_outerSize + outerChange; i++) m_innerNonZeros[i] = 0;
        }
        else if (innerChange < 0)
        {
            // Inner size decreased: allocate a new m_innerNonZeros
            m_innerNonZeros = static_cast<StorageIndex*>(std::malloc((m_outerSize + outerChange) * sizeof(StorageIndex)));
            if (!m_innerNonZeros)
                internal::throw_std_bad_alloc();
            for (Index i = 0; i < m_outerSize + (std::min)(outerChange, Index(0)); i++) m_innerNonZeros[i] = m_outerIndex[i + 1] - m_outerIndex[i];
            for (Index i = m_outerSize; i < m_outerSize + outerChange; i++) m_innerNonZeros[i] = 0;
        }

        // Change the m_innerNonZeros in case of a decrease of inner size
        if (m_innerNonZeros && innerChange < 0)
        {
            for (Index i = 0; i < m_outerSize + (std::min)(outerChange, Index(0)); i++)
            {
                StorageIndex& n = m_innerNonZeros[i];
                StorageIndex start = m_outerIndex[i];
                while (n > 0 && m_data.index(start + n - 1) >= newInnerSize) --n;
            }
        }

        m_innerSize = newInnerSize;

        // Re-allocate outer index structure if necessary
        if (outerChange == 0)
            return;

        StorageIndex* newOuterIndex = static_cast<StorageIndex*>(std::realloc(m_outerIndex, (m_outerSize + outerChange + 1) * sizeof(StorageIndex)));
        if (!newOuterIndex)
            internal::throw_std_bad_alloc();
        m_outerIndex = newOuterIndex;
        if (outerChange > 0)
        {
            StorageIndex lastIdx = m_outerSize == 0 ? 0 : m_outerIndex[m_outerSize];
            for (Index i = m_outerSize; i < m_outerSize + outerChange + 1; i++) m_outerIndex[i] = lastIdx;
        }
        m_outerSize += outerChange;
    }

    /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero.
      * 
      * This function does not free the currently allocated memory. To release as much as memory as possible,
      * call \code mat.data().squeeze(); \endcode after resizing it.
      * 
      * \sa reserve(), setZero()
      */
    void resize(Index rows, Index cols)
    {
        const Index outerSize = IsRowMajor ? rows : cols;
        m_innerSize = IsRowMajor ? cols : rows;
        m_data.clear();
        if (m_outerSize != outerSize || m_outerSize == 0)
        {
            std::free(m_outerIndex);
            m_outerIndex = static_cast<StorageIndex*>(std::malloc((outerSize + 1) * sizeof(StorageIndex)));
            if (!m_outerIndex)
                internal::throw_std_bad_alloc();

            m_outerSize = outerSize;
        }
        if (m_innerNonZeros)
        {
            std::free(m_innerNonZeros);
            m_innerNonZeros = 0;
        }
        memset(m_outerIndex, 0, (m_outerSize + 1) * sizeof(StorageIndex));
    }

    /** \internal
      * Resize the nonzero vector to \a size */
    void resizeNonZeros(Index size) { m_data.resize(size); }

    /** \returns a const expression of the diagonal coefficients. */
    const ConstDiagonalReturnType diagonal() const { return ConstDiagonalReturnType(*this); }

    /** \returns a read-write expression of the diagonal coefficients.
      * \warning If the diagonal entries are written, then all diagonal
      * entries \b must already exist, otherwise an assertion will be raised.
      */
    DiagonalReturnType diagonal() { return DiagonalReturnType(*this); }

    /** Default constructor yielding an empty \c 0 \c x \c 0 matrix */
    inline SparseMatrix() : m_outerSize(-1), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
    {
        check_template_parameters();
        resize(0, 0);
    }

    /** Constructs a \a rows \c x \a cols empty matrix */
    inline SparseMatrix(Index rows, Index cols) : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
    {
        check_template_parameters();
        resize(rows, cols);
    }

    /** Constructs a sparse matrix from the sparse expression \a other */
    template <typename OtherDerived>
    inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other) : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
    {
        EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
                            YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
        check_template_parameters();
        const bool needToTranspose = (Flags & RowMajorBit) != (internal::evaluator<OtherDerived>::Flags & RowMajorBit);
        if (needToTranspose)
            *this = other.derived();
        else
        {
#ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
            EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
#endif
            internal::call_assignment_no_alias(*this, other.derived());
        }
    }

    /** Constructs a sparse matrix from the sparse selfadjoint view \a other */
    template <typename OtherDerived, unsigned int UpLo>
    inline SparseMatrix(const SparseSelfAdjointView<OtherDerived, UpLo>& other) : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
    {
        check_template_parameters();
        Base::operator=(other);
    }

    /** Copy constructor (it performs a deep copy) */
    inline SparseMatrix(const SparseMatrix& other) : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
    {
        check_template_parameters();
        *this = other.derived();
    }

    /** \brief Copy constructor with in-place evaluation */
    template <typename OtherDerived>
    SparseMatrix(const ReturnByValue<OtherDerived>& other) : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
    {
        check_template_parameters();
        initAssignment(other);
        other.evalTo(*this);
    }

    /** \brief Copy constructor with in-place evaluation */
    template <typename OtherDerived>
    explicit SparseMatrix(const DiagonalBase<OtherDerived>& other) : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
    {
        check_template_parameters();
        *this = other.derived();
    }

    /** Swaps the content of two sparse matrices of the same type.
      * This is a fast operation that simply swaps the underlying pointers and parameters. */
    inline void swap(SparseMatrix& other)
    {
        //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
        std::swap(m_outerIndex, other.m_outerIndex);
        std::swap(m_innerSize, other.m_innerSize);
        std::swap(m_outerSize, other.m_outerSize);
        std::swap(m_innerNonZeros, other.m_innerNonZeros);
        m_data.swap(other.m_data);
    }

    /** Sets *this to the identity matrix.
      * This function also turns the matrix into compressed mode, and drop any reserved memory. */
    inline void setIdentity()
    {
        eigen_assert(rows() == cols() && "ONLY FOR SQUARED MATRICES");
        this->m_data.resize(rows());
        Eigen::Map<IndexVector>(this->m_data.indexPtr(), rows()).setLinSpaced(0, StorageIndex(rows() - 1));
        Eigen::Map<ScalarVector>(this->m_data.valuePtr(), rows()).setOnes();
        Eigen::Map<IndexVector>(this->m_outerIndex, rows() + 1).setLinSpaced(0, StorageIndex(rows()));
        std::free(m_innerNonZeros);
        m_innerNonZeros = 0;
    }
    inline SparseMatrix& operator=(const SparseMatrix& other)
    {
        if (other.isRValue())
        {
            swap(other.const_cast_derived());
        }
        else if (this != &other)
        {
#ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
            EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
#endif
            initAssignment(other);
            if (other.isCompressed())
            {
                internal::smart_copy(other.m_outerIndex, other.m_outerIndex + m_outerSize + 1, m_outerIndex);
                m_data = other.m_data;
            }
            else
            {
                Base::operator=(other);
            }
        }
        return *this;
    }

#ifndef EIGEN_PARSED_BY_DOXYGEN
    template <typename OtherDerived> inline SparseMatrix& operator=(const EigenBase<OtherDerived>& other) { return Base::operator=(other.derived()); }

    template <typename Lhs, typename Rhs> inline SparseMatrix& operator=(const Product<Lhs, Rhs, AliasFreeProduct>& other);
#endif  // EIGEN_PARSED_BY_DOXYGEN

    template <typename OtherDerived> EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other);

    friend std::ostream& operator<<(std::ostream& s, const SparseMatrix& m)
    {
        EIGEN_DBG_SPARSE(s << "Nonzero entries:\n"; if (m.isCompressed()) {
            for (Index i = 0; i < m.nonZeros(); ++i) s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
        } else {
            for (Index i = 0; i < m.outerSize(); ++i)
            {
                Index p = m.m_outerIndex[i];
                Index pe = m.m_outerIndex[i] + m.m_innerNonZeros[i];
                Index k = p;
                for (; k < pe; ++k) { s << "(" << m.m_data.value(k) << "," << m.m_data.index(k) << ") "; }
                for (; k < m.m_outerIndex[i + 1]; ++k) { s << "(_,_) "; }
            }
        } s << std::endl;
                         s << std::endl;
                         s << "Outer pointers:\n";
                         for (Index i = 0; i < m.outerSize(); ++i) { s << m.m_outerIndex[i] << " "; } s << " $" << std::endl;
                         if (!m.isCompressed()) {
                             s << "Inner non zeros:\n";
                             for (Index i = 0; i < m.outerSize(); ++i) { s << m.m_innerNonZeros[i] << " "; }
                             s << " $" << std::endl;
                         } s
                         << std::endl;);
        s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);
        return s;
    }

    /** Destructor */
    inline ~SparseMatrix()
    {
        std::free(m_outerIndex);
        std::free(m_innerNonZeros);
    }

    /** Overloaded for performance */
    Scalar sum() const;

#ifdef EIGEN_SPARSEMATRIX_PLUGIN
#include EIGEN_SPARSEMATRIX_PLUGIN
#endif

protected:
    template <typename Other> void initAssignment(const Other& other)
    {
        resize(other.rows(), other.cols());
        if (m_innerNonZeros)
        {
            std::free(m_innerNonZeros);
            m_innerNonZeros = 0;
        }
    }

    /** \internal
      * \sa insert(Index,Index) */
    EIGEN_DONT_INLINE Scalar& insertCompressed(Index row, Index col);

    /** \internal
      * A vector object that is equal to 0 everywhere but v at the position i */
    class SingletonVector
    {
        StorageIndex m_index;
        StorageIndex m_value;

    public:
        typedef StorageIndex value_type;
        SingletonVector(Index i, Index v) : m_index(convert_index(i)), m_value(convert_index(v)) {}

        StorageIndex operator[](Index i) const { return i == m_index ? m_value : 0; }
    };

    /** \internal
      * \sa insert(Index,Index) */
    EIGEN_DONT_INLINE Scalar& insertUncompressed(Index row, Index col);

public:
    /** \internal
      * \sa insert(Index,Index) */
    EIGEN_STRONG_INLINE Scalar& insertBackUncompressed(Index row, Index col)
    {
        const Index outer = IsRowMajor ? row : col;
        const Index inner = IsRowMajor ? col : row;

        eigen_assert(!isCompressed());
        eigen_assert(m_innerNonZeros[outer] <= (m_outerIndex[outer + 1] - m_outerIndex[outer]));

        Index p = m_outerIndex[outer] + m_innerNonZeros[outer]++;
        m_data.index(p) = convert_index(inner);
        return (m_data.value(p) = Scalar(0));
    }

protected:
    struct IndexPosPair
    {
        IndexPosPair(Index a_i, Index a_p) : i(a_i), p(a_p) {}
        Index i;
        Index p;
    };

    /** \internal assign \a diagXpr to the diagonal of \c *this
      * There are different strategies:
      *   1 - if *this is overwritten (Func==assign_op) or *this is empty, then we can work treat *this as a dense vector expression.
      *   2 - otherwise, for each diagonal coeff,
      *     2.a - if it already exists, then we update it,
      *     2.b - otherwise, if *this is uncompressed and that the current inner-vector has empty room for at least 1 element, then we perform an in-place insertion.
      *     2.c - otherwise, we'll have to reallocate and copy everything, so instead of doing so for each new element, it is recorded in a std::vector.
      *   3 - at the end, if some entries failed to be inserted in-place, then we alloc a new buffer, copy each chunk at the right position, and insert the new elements.
      * 
      * TODO: some piece of code could be isolated and reused for a general in-place update strategy.
      * TODO: if we start to defer the insertion of some elements (i.e., case 2.c executed once),
      *       then it *might* be better to disable case 2.b since they will have to be copied anyway.
      */
    template <typename DiagXpr, typename Func> void assignDiagonal(const DiagXpr diagXpr, const Func& assignFunc)
    {
        Index n = diagXpr.size();

        const bool overwrite = internal::is_same<Func, internal::assign_op<Scalar, Scalar>>::value;
        if (overwrite)
        {
            if ((this->rows() != n) || (this->cols() != n))
                this->resize(n, n);
        }

        if (m_data.size() == 0 || overwrite)
        {
            typedef Array<StorageIndex, Dynamic, 1> ArrayXI;
            this->makeCompressed();
            this->resizeNonZeros(n);
            Eigen::Map<ArrayXI>(this->innerIndexPtr(), n).setLinSpaced(0, StorageIndex(n) - 1);
            Eigen::Map<ArrayXI>(this->outerIndexPtr(), n + 1).setLinSpaced(0, StorageIndex(n));
            Eigen::Map<Array<Scalar, Dynamic, 1>> values = this->coeffs();
            values.setZero();
            internal::call_assignment_no_alias(values, diagXpr, assignFunc);
        }
        else
        {
            bool isComp = isCompressed();
            internal::evaluator<DiagXpr> diaEval(diagXpr);
            std::vector<IndexPosPair> newEntries;

            // 1 - try in-place update and record insertion failures
            for (Index i = 0; i < n; ++i)
            {
                internal::LowerBoundIndex lb = this->lower_bound(i, i);
                Index p = lb.value;
                if (lb.found)
                {
                    // the coeff already exists
                    assignFunc.assignCoeff(m_data.value(p), diaEval.coeff(i));
                }
                else if ((!isComp) && m_innerNonZeros[i] < (m_outerIndex[i + 1] - m_outerIndex[i]))
                {
                    // non compressed mode with local room for inserting one element
                    m_data.moveChunk(p, p + 1, m_outerIndex[i] + m_innerNonZeros[i] - p);
                    m_innerNonZeros[i]++;
                    m_data.value(p) = Scalar(0);
                    m_data.index(p) = StorageIndex(i);
                    assignFunc.assignCoeff(m_data.value(p), diaEval.coeff(i));
                }
                else
                {
                    // defer insertion
                    newEntries.push_back(IndexPosPair(i, p));
                }
            }
            // 2 - insert deferred entries
            Index n_entries = Index(newEntries.size());
            if (n_entries > 0)
            {
                Storage newData(m_data.size() + n_entries);
                Index prev_p = 0;
                Index prev_i = 0;
                for (Index k = 0; k < n_entries; ++k)
                {
                    Index i = newEntries[k].i;
                    Index p = newEntries[k].p;
                    internal::smart_copy(m_data.valuePtr() + prev_p, m_data.valuePtr() + p, newData.valuePtr() + prev_p + k);
                    internal::smart_copy(m_data.indexPtr() + prev_p, m_data.indexPtr() + p, newData.indexPtr() + prev_p + k);
                    for (Index j = prev_i; j < i; ++j) m_outerIndex[j + 1] += k;
                    if (!isComp)
                        m_innerNonZeros[i]++;
                    prev_p = p;
                    prev_i = i;
                    newData.value(p + k) = Scalar(0);
                    newData.index(p + k) = StorageIndex(i);
                    assignFunc.assignCoeff(newData.value(p + k), diaEval.coeff(i));
                }
                {
                    internal::smart_copy(m_data.valuePtr() + prev_p, m_data.valuePtr() + m_data.size(), newData.valuePtr() + prev_p + n_entries);
                    internal::smart_copy(m_data.indexPtr() + prev_p, m_data.indexPtr() + m_data.size(), newData.indexPtr() + prev_p + n_entries);
                    for (Index j = prev_i + 1; j <= m_outerSize; ++j) m_outerIndex[j] += n_entries;
                }
                m_data.swap(newData);
            }
        }
    }

private:
    static void check_template_parameters()
    {
        EIGEN_STATIC_ASSERT(NumTraits<StorageIndex>::IsSigned, THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE);
        EIGEN_STATIC_ASSERT((Options & (ColMajor | RowMajor)) == Options, INVALID_MATRIX_TEMPLATE_PARAMETERS);
    }

    struct default_prunning_func
    {
        default_prunning_func(const Scalar& ref, const RealScalar& eps) : reference(ref), epsilon(eps) {}
        inline bool operator()(const Index&, const Index&, const Scalar& value) const { return !internal::isMuchSmallerThan(value, reference, epsilon); }
        Scalar reference;
        RealScalar epsilon;
    };
};

namespace internal {

    template <typename InputIterator, typename SparseMatrixType, typename DupFunctor>
    void set_from_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat, DupFunctor dup_func)
    {
        enum
        {
            IsRowMajor = SparseMatrixType::IsRowMajor
        };
        typedef typename SparseMatrixType::Scalar Scalar;
        typedef typename SparseMatrixType::StorageIndex StorageIndex;
        SparseMatrix<Scalar, IsRowMajor ? ColMajor : RowMajor, StorageIndex> trMat(mat.rows(), mat.cols());

        if (begin != end)
        {
            // pass 1: count the nnz per inner-vector
            typename SparseMatrixType::IndexVector wi(trMat.outerSize());
            wi.setZero();
            for (InputIterator it(begin); it != end; ++it)
            {
                eigen_assert(it->row() >= 0 && it->row() < mat.rows() && it->col() >= 0 && it->col() < mat.cols());
                wi(IsRowMajor ? it->col() : it->row())++;
            }

            // pass 2: insert all the elements into trMat
            trMat.reserve(wi);
            for (InputIterator it(begin); it != end; ++it) trMat.insertBackUncompressed(it->row(), it->col()) = it->value();

            // pass 3:
            trMat.collapseDuplicates(dup_func);
        }

        // pass 4: transposed copy -> implicit sorting
        mat = trMat;
    }

}  // namespace internal

/** Fill the matrix \c *this with the list of \em triplets defined by the iterator range \a begin - \a end.
  *
  * A \em triplet is a tuple (i,j,value) defining a non-zero element.
  * The input list of triplets does not have to be sorted, and can contains duplicated elements.
  * In any case, the result is a \b sorted and \b compressed sparse matrix where the duplicates have been summed up.
  * This is a \em O(n) operation, with \em n the number of triplet elements.
  * The initial contents of \c *this is destroyed.
  * The matrix \c *this must be properly resized beforehand using the SparseMatrix(Index,Index) constructor,
  * or the resize(Index,Index) method. The sizes are not extracted from the triplet list.
  *
  * The \a InputIterators value_type must provide the following interface:
  * \code
  * Scalar value() const; // the value
  * Scalar row() const;   // the row index i
  * Scalar col() const;   // the column index j
  * \endcode
  * See for instance the Eigen::Triplet template class.
  *
  * Here is a typical usage example:
  * \code
    typedef Triplet<double> T;
    std::vector<T> tripletList;
    tripletList.reserve(estimation_of_entries);
    for(...)
    {
      // ...
      tripletList.push_back(T(i,j,v_ij));
    }
    SparseMatrixType m(rows,cols);
    m.setFromTriplets(tripletList.begin(), tripletList.end());
    // m is ready to go!
  * \endcode
  *
  * \warning The list of triplets is read multiple times (at least twice). Therefore, it is not recommended to define
  * an abstract iterator over a complex data-structure that would be expensive to evaluate. The triplets should rather
  * be explicitly stored into a std::vector for instance.
  */
template <typename Scalar, int _Options, typename _StorageIndex>
template <typename InputIterators>
void SparseMatrix<Scalar, _Options, _StorageIndex>::setFromTriplets(const InputIterators& begin, const InputIterators& end)
{
    internal::set_from_triplets<InputIterators, SparseMatrix<Scalar, _Options, _StorageIndex>>(begin, end, *this, internal::scalar_sum_op<Scalar, Scalar>());
}

/** The same as setFromTriplets but when duplicates are met the functor \a dup_func is applied:
  * \code
  * value = dup_func(OldValue, NewValue)
  * \endcode 
  * Here is a C++11 example keeping the latest entry only:
  * \code
  * mat.setFromTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; });
  * \endcode
  */
template <typename Scalar, int _Options, typename _StorageIndex>
template <typename InputIterators, typename DupFunctor>
void SparseMatrix<Scalar, _Options, _StorageIndex>::setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func)
{
    internal::set_from_triplets<InputIterators, SparseMatrix<Scalar, _Options, _StorageIndex>, DupFunctor>(begin, end, *this, dup_func);
}

/** \internal */
template <typename Scalar, int _Options, typename _StorageIndex>
template <typename DupFunctor>
void SparseMatrix<Scalar, _Options, _StorageIndex>::collapseDuplicates(DupFunctor dup_func)
{
    eigen_assert(!isCompressed());
    // TODO, in practice we should be able to use m_innerNonZeros for that task
    IndexVector wi(innerSize());
    wi.fill(-1);
    StorageIndex count = 0;
    // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers
    for (Index j = 0; j < outerSize(); ++j)
    {
        StorageIndex start = count;
        Index oldEnd = m_outerIndex[j] + m_innerNonZeros[j];
        for (Index k = m_outerIndex[j]; k < oldEnd; ++k)
        {
            Index i = m_data.index(k);
            if (wi(i) >= start)
            {
                // we already meet this entry => accumulate it
                m_data.value(wi(i)) = dup_func(m_data.value(wi(i)), m_data.value(k));
            }
            else
            {
                m_data.value(count) = m_data.value(k);
                m_data.index(count) = m_data.index(k);
                wi(i) = count;
                ++count;
            }
        }
        m_outerIndex[j] = start;
    }
    m_outerIndex[m_outerSize] = count;

    // turn the matrix into compressed form
    std::free(m_innerNonZeros);
    m_innerNonZeros = 0;
    m_data.resize(m_outerIndex[m_outerSize]);
}

template <typename Scalar, int _Options, typename _StorageIndex>
template <typename OtherDerived>
EIGEN_DONT_INLINE SparseMatrix<Scalar, _Options, _StorageIndex>& SparseMatrix<Scalar, _Options, _StorageIndex>::
operator=(const SparseMatrixBase<OtherDerived>& other)
{
    EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
                        YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)

#ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
    EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
#endif

    const bool needToTranspose = (Flags & RowMajorBit) != (internal::evaluator<OtherDerived>::Flags & RowMajorBit);
    if (needToTranspose)
    {
#ifdef EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN
        EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN
#endif
        // two passes algorithm:
        //  1 - compute the number of coeffs per dest inner vector
        //  2 - do the actual copy/eval
        // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed
        typedef typename internal::nested_eval<OtherDerived, 2, typename internal::plain_matrix_type<OtherDerived>::type>::type OtherCopy;
        typedef typename internal::remove_all<OtherCopy>::type _OtherCopy;
        typedef internal::evaluator<_OtherCopy> OtherCopyEval;
        OtherCopy otherCopy(other.derived());
        OtherCopyEval otherCopyEval(otherCopy);

        SparseMatrix dest(other.rows(), other.cols());
        Eigen::Map<IndexVector>(dest.m_outerIndex, dest.outerSize()).setZero();

        // pass 1
        // FIXME the above copy could be merged with that pass
        for (Index j = 0; j < otherCopy.outerSize(); ++j)
            for (typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it) ++dest.m_outerIndex[it.index()];

        // prefix sum
        StorageIndex count = 0;
        IndexVector positions(dest.outerSize());
        for (Index j = 0; j < dest.outerSize(); ++j)
        {
            StorageIndex tmp = dest.m_outerIndex[j];
            dest.m_outerIndex[j] = count;
            positions[j] = count;
            count += tmp;
        }
        dest.m_outerIndex[dest.outerSize()] = count;
        // alloc
        dest.m_data.resize(count);
        // pass 2
        for (StorageIndex j = 0; j < otherCopy.outerSize(); ++j)
        {
            for (typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it)
            {
                Index pos = positions[it.index()]++;
                dest.m_data.index(pos) = j;
                dest.m_data.value(pos) = it.value();
            }
        }
        this->swap(dest);
        return *this;
    }
    else
    {
        if (other.isRValue())
        {
            initAssignment(other.derived());
        }
        // there is no special optimization
        return Base::operator=(other.derived());
    }
}

template <typename _Scalar, int _Options, typename _StorageIndex>
typename SparseMatrix<_Scalar, _Options, _StorageIndex>::Scalar& SparseMatrix<_Scalar, _Options, _StorageIndex>::insert(Index row, Index col)
{
    eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols());

    const Index outer = IsRowMajor ? row : col;
    const Index inner = IsRowMajor ? col : row;

    if (isCompressed())
    {
        if (nonZeros() == 0)
        {
            // reserve space if not already done
            if (m_data.allocatedSize() == 0)
                m_data.reserve(2 * m_innerSize);

            // turn the matrix into non-compressed mode
            m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
            if (!m_innerNonZeros)
                internal::throw_std_bad_alloc();

            memset(m_innerNonZeros, 0, (m_outerSize) * sizeof(StorageIndex));

            // pack all inner-vectors to the end of the pre-allocated space
            // and allocate the entire free-space to the first inner-vector
            StorageIndex end = convert_index(m_data.allocatedSize());
            for (Index j = 1; j <= m_outerSize; ++j) m_outerIndex[j] = end;
        }
        else
        {
            // turn the matrix into non-compressed mode
            m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
            if (!m_innerNonZeros)
                internal::throw_std_bad_alloc();
            for (Index j = 0; j < m_outerSize; ++j) m_innerNonZeros[j] = m_outerIndex[j + 1] - m_outerIndex[j];
        }
    }

    // check whether we can do a fast "push back" insertion
    Index data_end = m_data.allocatedSize();

    // First case: we are filling a new inner vector which is packed at the end.
    // We assume that all remaining inner-vectors are also empty and packed to the end.
    if (m_outerIndex[outer] == data_end)
    {
        eigen_internal_assert(m_innerNonZeros[outer] == 0);

        // pack previous empty inner-vectors to end of the used-space
        // and allocate the entire free-space to the current inner-vector.
        StorageIndex p = convert_index(m_data.size());
        Index j = outer;
        while (j >= 0 && m_innerNonZeros[j] == 0) m_outerIndex[j--] = p;

        // push back the new element
        ++m_innerNonZeros[outer];
        m_data.append(Scalar(0), inner);

        // check for reallocation
        if (data_end != m_data.allocatedSize())
        {
            // m_data has been reallocated
            //  -> move remaining inner-vectors back to the end of the free-space
            //     so that the entire free-space is allocated to the current inner-vector.
            eigen_internal_assert(data_end < m_data.allocatedSize());
            StorageIndex new_end = convert_index(m_data.allocatedSize());
            for (Index k = outer + 1; k <= m_outerSize; ++k)
                if (m_outerIndex[k] == data_end)
                    m_outerIndex[k] = new_end;
        }
        return m_data.value(p);
    }

    // Second case: the next inner-vector is packed to the end
    // and the current inner-vector end match the used-space.
    if (m_outerIndex[outer + 1] == data_end && m_outerIndex[outer] + m_innerNonZeros[outer] == m_data.size())
    {
        eigen_internal_assert(outer + 1 == m_outerSize || m_innerNonZeros[outer + 1] == 0);

        // add space for the new element
        ++m_innerNonZeros[outer];
        m_data.resize(m_data.size() + 1);

        // check for reallocation
        if (data_end != m_data.allocatedSize())
        {
            // m_data has been reallocated
            //  -> move remaining inner-vectors back to the end of the free-space
            //     so that the entire free-space is allocated to the current inner-vector.
            eigen_internal_assert(data_end < m_data.allocatedSize());
            StorageIndex new_end = convert_index(m_data.allocatedSize());
            for (Index k = outer + 1; k <= m_outerSize; ++k)
                if (m_outerIndex[k] == data_end)
                    m_outerIndex[k] = new_end;
        }

        // and insert it at the right position (sorted insertion)
        Index startId = m_outerIndex[outer];
        Index p = m_outerIndex[outer] + m_innerNonZeros[outer] - 1;
        while ((p > startId) && (m_data.index(p - 1) > inner))
        {
            m_data.index(p) = m_data.index(p - 1);
            m_data.value(p) = m_data.value(p - 1);
            --p;
        }

        m_data.index(p) = convert_index(inner);
        return (m_data.value(p) = Scalar(0));
    }

    if (m_data.size() != m_data.allocatedSize())
    {
        // make sure the matrix is compatible to random un-compressed insertion:
        m_data.resize(m_data.allocatedSize());
        this->reserveInnerVectors(Array<StorageIndex, Dynamic, 1>::Constant(m_outerSize, 2));
    }

    return insertUncompressed(row, col);
}

template <typename _Scalar, int _Options, typename _StorageIndex>
EIGEN_DONT_INLINE typename SparseMatrix<_Scalar, _Options, _StorageIndex>::Scalar& SparseMatrix<_Scalar, _Options, _StorageIndex>::insertUncompressed(Index row,
                                                                                                                                                      Index col)
{
    eigen_assert(!isCompressed());

    const Index outer = IsRowMajor ? row : col;
    const StorageIndex inner = convert_index(IsRowMajor ? col : row);

    Index room = m_outerIndex[outer + 1] - m_outerIndex[outer];
    StorageIndex innerNNZ = m_innerNonZeros[outer];
    if (innerNNZ >= room)
    {
        // this inner vector is full, we need to reallocate the whole buffer :(
        reserve(SingletonVector(outer, std::max<StorageIndex>(2, innerNNZ)));
    }

    Index startId = m_outerIndex[outer];
    Index p = startId + m_innerNonZeros[outer];
    while ((p > startId) && (m_data.index(p - 1) > inner))
    {
        m_data.index(p) = m_data.index(p - 1);
        m_data.value(p) = m_data.value(p - 1);
        --p;
    }
    eigen_assert((p <= startId || m_data.index(p - 1) != inner) && "you cannot insert an element that already exists, you must call coeffRef to this end");

    m_innerNonZeros[outer]++;

    m_data.index(p) = inner;
    return (m_data.value(p) = Scalar(0));
}

template <typename _Scalar, int _Options, typename _StorageIndex>
EIGEN_DONT_INLINE typename SparseMatrix<_Scalar, _Options, _StorageIndex>::Scalar& SparseMatrix<_Scalar, _Options, _StorageIndex>::insertCompressed(Index row,
                                                                                                                                                    Index col)
{
    eigen_assert(isCompressed());

    const Index outer = IsRowMajor ? row : col;
    const Index inner = IsRowMajor ? col : row;

    Index previousOuter = outer;
    if (m_outerIndex[outer + 1] == 0)
    {
        // we start a new inner vector
        while (previousOuter >= 0 && m_outerIndex[previousOuter] == 0)
        {
            m_outerIndex[previousOuter] = convert_index(m_data.size());
            --previousOuter;
        }
        m_outerIndex[outer + 1] = m_outerIndex[outer];
    }

    // here we have to handle the tricky case where the outerIndex array
    // starts with: [ 0 0 0 0 0 1 ...] and we are inserted in, e.g.,
    // the 2nd inner vector...
    bool isLastVec = (!(previousOuter == -1 && m_data.size() != 0)) && (std::size_t(m_outerIndex[outer + 1]) == m_data.size());

    std::size_t startId = m_outerIndex[outer];
    // FIXME let's make sure sizeof(long int) == sizeof(std::size_t)
    std::size_t p = m_outerIndex[outer + 1];
    ++m_outerIndex[outer + 1];

    double reallocRatio = 1;
    if (m_data.allocatedSize() <= m_data.size())
    {
        // if there is no preallocated memory, let's reserve a minimum of 32 elements
        if (m_data.size() == 0)
        {
            m_data.reserve(32);
        }
        else
        {
            // we need to reallocate the data, to reduce multiple reallocations
            // we use a smart resize algorithm based on the current filling ratio
            // in addition, we use double to avoid integers overflows
            double nnzEstimate = double(m_outerIndex[outer]) * double(m_outerSize) / double(outer + 1);
            reallocRatio = (nnzEstimate - double(m_data.size())) / double(m_data.size());
            // furthermore we bound the realloc ratio to:
            //   1) reduce multiple minor realloc when the matrix is almost filled
            //   2) avoid to allocate too much memory when the matrix is almost empty
            reallocRatio = (std::min)((std::max)(reallocRatio, 1.5), 8.);
        }
    }
    m_data.resize(m_data.size() + 1, reallocRatio);

    if (!isLastVec)
    {
        if (previousOuter == -1)
        {
            // oops wrong guess.
            // let's correct the outer offsets
            for (Index k = 0; k <= (outer + 1); ++k) m_outerIndex[k] = 0;
            Index k = outer + 1;
            while (m_outerIndex[k] == 0) m_outerIndex[k++] = 1;
            while (k <= m_outerSize && m_outerIndex[k] != 0) m_outerIndex[k++]++;
            p = 0;
            --k;
            k = m_outerIndex[k] - 1;
            while (k > 0)
            {
                m_data.index(k) = m_data.index(k - 1);
                m_data.value(k) = m_data.value(k - 1);
                k--;
            }
        }
        else
        {
            // we are not inserting into the last inner vec
            // update outer indices:
            Index j = outer + 2;
            while (j <= m_outerSize && m_outerIndex[j] != 0) m_outerIndex[j++]++;
            --j;
            // shift data of last vecs:
            Index k = m_outerIndex[j] - 1;
            while (k >= Index(p))
            {
                m_data.index(k) = m_data.index(k - 1);
                m_data.value(k) = m_data.value(k - 1);
                k--;
            }
        }
    }

    while ((p > startId) && (m_data.index(p - 1) > inner))
    {
        m_data.index(p) = m_data.index(p - 1);
        m_data.value(p) = m_data.value(p - 1);
        --p;
    }

    m_data.index(p) = inner;
    return (m_data.value(p) = Scalar(0));
}

namespace internal {

    template <typename _Scalar, int _Options, typename _StorageIndex>
    struct evaluator<SparseMatrix<_Scalar, _Options, _StorageIndex>> : evaluator<SparseCompressedBase<SparseMatrix<_Scalar, _Options, _StorageIndex>>>
    {
        typedef evaluator<SparseCompressedBase<SparseMatrix<_Scalar, _Options, _StorageIndex>>> Base;
        typedef SparseMatrix<_Scalar, _Options, _StorageIndex> SparseMatrixType;
        evaluator() : Base() {}
        explicit evaluator(const SparseMatrixType& mat) : Base(mat) {}
    };

}  // namespace internal

}  // end namespace Eigen

#endif  // EIGEN_SPARSEMATRIX_H
